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Abstract 

Recently launched high precision Synthetic Aperture Radar (SAR) satellites such as TerraSAR- 
X, COSMO-SkyMed, etc. present a high potential for better observation and characterization of the 
cryosphere. This study introduces a new approach using high frequency (X-band) SAR data and an 
Electromagnetic Backscattering Model (EBM) to constrain the detailed snowpack model Crocus. A 
snowpack EBM based on radiative transfer theory, previously used for C-band applications, is adapted 
for the X-band. From measured or simulated snowpack stratigraphic profiles consisting of snow optical 
grain radius and density, this forward model calculates the backscattering coefficient a for different 
polarimetric channels. The output result is then compared with spaceborne TerraSAR-X acquisitions 
to evaluate the forward model. Next, from the EBM, the adjoint operator is developed and used in 
a variational analysis scheme in order to minimize the discrepancies between simulations and SAR 
observations. A time series of TerraSAR-X acquisitions and in-situ measurements on the Argentiere 
glacier (Mont-Blanc massif, French Alps) are used to evaluate the EBM and the data assimilation 
scheme. Results indicate that snow stratigraphic profiles obtained after the analysis process show a 
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closer agreement with the measured ones than the initial ones, and therefore demonstrate the high 
potential of assimilating SAR data to model of snow evolution. 

Index Terms 

Remote sensing, electromagnetic backscattering model, snow grain size, snow density, radar (SAR), 
data analysis. 

I. Introduction 

Snowpack characterization has become a critical issue in the present context of climate change. 
Estimating some of the properties of a snowpack, like its density and grain size distribution 
will provide great benefit to snow forecasting, prevision of natural hazard, like snow avalanche 
warning, and economic arrangements related to tourism and winter sports. Due to its imaging 
capabilities over large areas, unaffected by weather and day-night conditions, Synthetic Aperture 
Radar (SAR) is an important tool for snowpack characterization in a natural environment. 
Moreover, the high penetration depth of radar electromagnetic waves allow us to retrieve the 
information inside the volume of the snowpack. Over the past decade, the large availability of 
L and C-band SAR data provided by various spaceborne sensors, like ALOS PALSAR, ERS-1, 
ENVISAT, led to many studies on the characterization of snowpack properties [[TJ, [|2j. 

A new generation of X-band (8-12GHz) SAR systems, and in the near future Ku-band (12- 
18GHz), with high image resolution, short revisit time will provide improved information that 
might be used to characterize and monitor snowpack. In this context, it is necessary to develop 
a compatible EBM accounting for electromagnetic waves (EMW) propagation and scattering at 
high frequencies (X and Ku-bands) through a multilayer snowpack. Some backscattering models 
at L and C-band frequencies have been introduced in [|2j, (3j. These models simulate the loss 
of EMW energy while propagating through dense media by solving the Radiative Transfer (RT) 
differential equation HJ. In order to introduce coherent recombination effects in the RT coherent 
model, Wang, et al. [|5) applied the Strong Fluctuation Theory (SFT) introduced by Stogryn (6j to 
calculate the effective permittivity of each snow layer, in which the correlation among particles 
was taken into account. The scattering and absorption mechanisms in the EBM are simulated 
using the Rayleigh scattering model due to the snow grain size being in this study is much 
smaller than the carrier wavelength. 
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In this paper, the snowpack backscattering model initially developed in [|2| is adapted for 
X-band and higher frequencies, in the case of a dry snow medium. The adaptation consists of 
updating the IEM introduced by Fung et. al. in 1992 fTJ by a newer version published in 2004 [8|, 
which allows the calculation of surface and ground backscattering components for X-band and 
higher frequencies. Meanwhile the modeling of volume backscattering of the existed model, 
which is based on solving the Vector Radiative Transfer equation and Rayleigh scattering model, 
is compatible for X and Ku-bands. From the physical features of each snow layer (optical grain 
radius, density, thickness) and for given SAR acquisition conditions (frequency, incidence angle), 
the model calculates the total backscattering coefficient cr° for different polarization channels 
and their vertical distribution within the snowpack. Next, the snowpack profiles generated by 
the detailed snowpack model Crocus using downscaled meteorological fields from the SAFRAN 
analysis [|9)-[ 1 1 1, are constrained using the SAR image data and EBM simulations. In this study, 
the number of observable, i.e. the SAR backscattering coefficients, being much smaller than the 
number of unknown parameters, i.e. the snow cover properties, a classical estimation approach 
based on the use of an inverse problem would reveal totally inefficient. Instead, an adjoint operator 
of the direct EBM is developed to be used in a assimilation scheme. A variational assimilation 
method allows the integration of the observation data into a set of initial guess parameters 
through a direct model, and therefore can constrain these parameters without explicitly inverting 



the model. In our study, the three-dimensional variational analysis (3D-VAR) method [12| is 
implemented. Finally, a time series of TerraSAR-X acquisitions on the mountainous region of 
the French-Alps is used to evaluate the model and the data assimilation process. The Argentiere 
glacier area has been chosen for the case study due to its large, uniformly snow-covered surface 
area. Some in-situ measurements on this area are also available at the same timeline of SAR 
acquisitions and therefore are used to evaluate the EBM and the performance of the data 
assimilation scheme. 

Details of the EBM equations and its the physical and mathematical hypothesis are presented 
in section II. An introduction to the Crocus detailed snowpack model and to the detailed imple- 
mentation of the 3D-VAR scheme are described in section III. Section IV shows a description 
of TerraSAR-X acquisition parameters, as well as results of the case study on the Argentiere 
glacier. 
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Fig. 1. Main backscattering mechanisms occurring within a multilayer snowpack that can be simulated using the RT theory at 
order 1: air-snow reflection (M as ), volume scattering (M„ i) and reflection over the ground (M 9 ). 



II. EMW BACKSCATTERING MODEL 

A. Main components of the total backscattering coefficient 

The Stoke vector, which contains the incoherent information related to the polarization of an 
electromagnetic wave, can be expressed as: 
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electric field [13|, and (.) represents the expectation operator. 

For given acquisition conditions, the Stoke vector scattered by a medium, g s , can be related 
to the incident one, g i5 by a Mueller matrix M as g s = Mgj with: 
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where M u = cr° v and M 22 = a® h represent the co-polarized backscattering coefficients and 
M 33 = Xe(pthh) and M 34 = $srn(a® vhh ) are correlation terms. Due to the reflection symmetry, 
the cross-polarization coefficients of the matrix a pppq are equal to zero p"3|, and the rest of the 
elements of M are null. 

The solution of the RT equation at order 1 provides a total backscattered information from 
a snowpack that consists of a combination of five scattering mechanisms: reflection at the 
surface air-snow interface, volume scattering, volume-ground and ground-volume interactions, 
and reflection over the ground | fl4| . Due to their small amplitude, the volume-ground and ground- 
volume contributions can be neglected [15]. The illustration of the three other mechanisms is 
shown in Fig.[T] The expression of the total polarimetric backscattered information can be written 
using the Mueller matrix corresponding to each mechanism: 

M snow = M as + M vol + M g (3) 

The surface and ground backscattering are modeled using the IEM introduced by Fung 
et. al. ][8j, whereas the volume contribution is calculated using the Vector Radiative Transfer 
equation. 

B. Surface backscattering 

The matrix M as represents the second order polarimetric response backscattered by the air- 
snow interface. Its elements can be calculated from the surface roughness parameters, e.g. its 
correlation function w(x) and its root mean square (rms) height cr^, the incidence angle 9q and 
the emitted EM wave frequency / using the Integration Equation Model (IEM) pj. According 
to the IEM, the reflectivity may be expressed as: 

a° pq = ^ exp(-2A; 2 ^ cos 2 B ) £ 1 °- } ' (4) 

where p and q are equal to h or v, indicating a horizontal or vertical polarization, k = 
represents the wave number. The detailed mathematical expressions of the surface spectrum W n 
and the Fresnel reflection/transmission factor \I™ q \ can be found in ][8j. 

C. Volume backscattering 

The volume backscattering ~M vo i is deduced from the loss of EMW intensity during propagation 
through a multilayer snowpack, which can be categorized into 4 types: related to transmission 
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between two layers, absorption by the snow particles, scattering and coherent recombination. 
The amplitude of each mechanism depends largely on the dielectric properties of the snowpack 
medium. Therefore the permittivity of each layer, which characterizes its dielectric properties, 
needs to be calculated first: 

1) Dry snow permittivity: Dry snow is considered as a dense and heterogeneous medium 
with strong variations of various physical properties such as grain size, density, thickness. 
Therefore the variance of permittivity across a snow layer is relatively high. Several snowpack 
characterization methods JT6| are largely based on the assumption that the scattering losses due 
to the correlation of EMW are negligible. However, at high frequency, the snowpack structure 
becomes bigger compared to the wavelength of X or Ku-band EM waves (T7J. The correlation 
between particles can no longer be ignored. The Strong Fluctuation Theory (SFT) introduced 
by Stogryn j6j can model the permittivity of such medium by using the effective permittivity 
(e e ff) that takes into account the scattering effect among ice particles at high frequencies. The 
expression of e e ff using SFT is as follows [jTVj : 



z eff 



(5) 



where e g and 8 €g are the quasi-static permittivity and its variance, k the wave number and L the 
correlation length, which is proportional to the average snow grain size and the snow density of 
the medium. 

2) Transmission between two layers: The snowpack consists of layers with different physical 
properties. Therefore the model needs to take into account the energy loss due to transmission 
between two layers. With the assumption of a smooth interface between two layers, the Fresnel 
transmission can be used. It is expressed through a matrix as follows fl4j: 
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represents the Fresnel transmission coefficients of pp channel, whereas gk(k-i) 
and h k (k-i) are the terms of Mueller matrix related to the co-polarized correlation |2j: 

COS gfc-l yw.tw ,hh* \ „_j t COS Afc-l cv, /.w .Wi* \ / 7 n 



3 J 77ze attenuation: The particles in a snowpack are generally considered as spheres [(3), [ 15 1, 



p8| . Due to the spherical symmetry of the particle shape, the extinction of a wave propagating 
through the snowpack is independent of the polarization and may hence be represented by a 
scalar coefficient. The extinction is composed of an absorption and a scattering terms: 

K e = K a + k s (8) 



It can also be computed through the effective permitivity e e // |17|: 

n e = 2k Im y/e^ff (9) 

The attenuation matrix represents the gradual loss in EMW intensity while penetrating through 
a multilayer snowpack, composed of layers with different physical properties. It takes into account 
the energy loss by absorption and scattering mechanisms based on the extinction coefficient K e 
and thickness d of the layer, as well as the loss by transmission effect while an EM propagate 
through different layers: 

k f n i d i \ 

Mt down (k) = J|exp e — Tj(j_!) (10) 

i=i V costly 

Mt up {k) = nT ( i-i)i exp (-^r) (11) 

The Attdown is the intensity loss when propagating from the surface to layer k, whereas Att up 
represents the intensity loss from layer k to the surface. The exponential factor, which takes 
into account the gradual loss of energy throughout the layer, is deduced from the basic radiative 
transfer equation dl = In e dr where r = djcosQ. 

4) Scattering by the particles: The phase matrix P fc under the hypothesis of spherical particles 
has the form shown in ([2]) where the cross-polarization terms P 12 and P 2 i are nuil - I n tne 
backscattering case, with the assumption of spherical particles, the SFT phase matrix can be 
simplified to P k = |^/ 4 where J 4 is the (4x4) identity matrix [19|. 
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5) Calculation of the volume backscattering: If we consider a snowpack made of n distinct 
layers, where 0& is the incidence angle and d k is the thickness of layer k, the total contribution 
of the volume backscattering mechanism M ra ; can be written as follows: 



M voi = 47rcos6»o Att up (/c - l)T (fe _i) fe 
fc=i 

2K*d k \ 
cos 9 k ) -nk 



1 — exp 



2k* 



P*T fc(fc _ 1) Att li0Wl (A; - 1) (12) 



Z). Ground backscattering 

The backscattering M g of the snow-ground interface may be computed as: 

M g = cos6 Att up (n)^^-AU down (n) (13) 

cos 6 n 

where R(0 n ) represents the contribution of the underlying ground surface backscattering and can 
be determined using the IEM. 

III. 3D- VAR DATA ASSIMILATION 

A. The detailed snowpack model Crocus 

Crocus is a one-dimensional numerical model simulating the thermodynamic balance of energy 
and mass of the snowpack. Its main objective is to describe in detail the evolution of internal 
snowpack properties based on the description of the evolution of morphological properties of 
snow grains during their metamorphism. Fig. [2] describes the general scheme of Crocus. It takes 
as input the meteorological variables air temperature, relative air humidity, wind speed, solar 
radiation, long wave radiation, amount and phase of precipitation. When it is used in the French 
mountain ranges (Alps, Pyrenees and Corsica), these quantities are commonly provided by the 
SAFRAN system, which combines ground-based and radiosondes observations with an a priori 
estimate of meteorological conditions from a numerical weather prediction (NWP) model (5), 



[ 1 1 . The output includes the scalar physical properties of the snowpack (snow depth, snow 
water equivalent, surface temperature, albedo, . . . ) along with the internal physical properties 
for each layer (density, thickness, optical grain radius, . . . ). SAFRAN meteorological fields are 
assumed to be homogeneous within a given mountrain range and provide a description of the 
altitude dependency of meteorological variables by steps of 300 m elevation ||9), [10|. 
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Fig. 2. General scheme of Crocus processes and variables 



Here we use the latest version of the detailed snowpack model Crocus, recently incorporated 



in the land surface scheme ISBA within the SURFEX interface [11 1. Among other advantages 
over previous versions of Crocus, this allows seamless coupling of the snowpack to the state of 
the underlying ground. 

B. Method introduction 

Variational assimilation aims to integrate observation data into guess parameters through the 
use of an observation operator. It is widely used in meteorological studies in order to relate 



observations, measurements and modeling aspects [20|. The method concentrates on searching a 
solution that minimizes simultaneously the distance between observations and simulation results 
and the distance between initial guess variables and the analysed variables. A scheme of this 
process is presented in Fig. [3] In this part of the paper, the output of our EBM in the previous 
section, such as backscattering coefficient of HH and VV polarizations, are used as elements of 
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I 

Snow backscattering 
model H(x) = o slm 




Fig. 3. Global scheme of the data assimilation used in this study. The input of the process are the SAR reflectivities, a , 
(observation) and the snowpack stratigraphic profile calculated by Crocus (guess). The output is the assimilated snowpack profile 
x that minimizes the cost function. 

the observation operator H ebm (x): 

H efem (x) = vec(M snow ) (14) 

where x represents the set of variables required to describe the snowpack properties. 

The 3D-VAR Jl2| algorithm is based on the minimization of a cost function J(x), defined as: 

J(x) = (x - x^B-^x - x g ) + (y obs - H eWt (x))*R- 1 (y - H e6m (x)) (15) 

where x is called the state vector, and can be modified after each iteration of minimization, x 5 
is the initial guess of the state vector and remains constant during the whole process. Therefore 
|| x — x g || 2 serves as a distance between the modified profile and the starting point. The observed 
polarimetric response, y obs , is denoted similarly to the calibrated values of the backscattering 
coefficients a . Therefore, ||y — H e f, m (x)|| 2 represents the distance between simulated and ob- 
served quantities in the observation space. The process also requires the estimation of the error 
covariance matrices of observations/simulations R and of the model B, the guess error variance. 
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C. Adjoint operator and minimization algorithm 

In order to minimize the cost function J, one needs to calculate its gradient: 

V J(x) = ^ = 2B- 1 (x - x 9 ) - 2VH* 6m (x)R- 1 (y obs - H e6m (x)) (16) 

If the model is denoted H efem : B — > 1Z, with B and 1Z are the domain of definition of x 
and y, then the function VH* 6m satisfying: Vx,y, (VH* ebm y,x) B = (y, VH efem x)7e is the adjoint 
operator of H ebm . 

Once the adjoint operator is developed, the minimization of J can be achieved using a gradient 
descent algorithm. Each iteration consists of modifying the vector x by a factor according to the 
Newton method: 

x„+i = x n + (V 2 J(x„))" 1 VJ(x n ) (17) 
where V 2 J(x n ) is the gradient of second order (hessian) of J: 

V 2 J = 2B- 1 + 2VH* fem R- 1 VH efm (18) 

D. Discussion on the assimilation method 

In general, the aim of modeling the relation between the elements of natural environment and 
the observations measured by special equipments (such as SAR or optical sensors) is to try to 
inverse the model and estimate the variables of environment using the observations. However, 
such problems often lead to the need to resolve a underdetermined system, which means the 
number of unknown is higher than the number of equations. 

In our case, the length of the input state vector x can reach 100 (in the case of snowpack 
with 50 layers, which is frequently generated by Crocus), meanwhile the output of the model 
consists of only the backscattering coefficients corresponding to polarimetric channels of SAR 
data. Therefore the realization of an inverse model is theoretically impossible. 

The data analysis method, on the other hand, requires a vector of guess variables relatively 
close to the actual values, allowing to add an a priori information. The snowpack variables 
calculated by Crocus are used as guess in our assimilation scheme. The fundamental goal is to 
try to modify the initial guess variables, based on balancing the errors of guess, modeling and 
measurements. It should be noted that the problem stays underdetermined, the analysis scheme 
only serves as a method to improve the initial guess variables using the new observations from 
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TABLE I 

TerraSAR-X acquisitions parameters 



Parameter 


Value 


TerraSAR-X products 


Single Look Complex Image 


Frequency (GHz) 


9.65 


Channels 


HH 


Incidence angle (deg) 


37.9892 


Mode 


Descending 


Acquisition dates 


6 Jan 2009, 17 Jan 2009, 28 Jan 2009 
8 Feb 2009, 19 Feb 2009, 2 Mar 2009 
13 Mar 2009, 24 Mar 2009 


Range resolution (m) 


1.477 


Azimuth resolution (m) 


2.44 


Calibration gain (dB) 


49.6802 



SAR data. The quality of improvement is based on the estimation of the initial guess vector x g 
and the precision of the EBM. 

IV. Case Study: Argentiere Glacier 

A. Data 

For this study, a time series of TerraSAR-X descending acquisitions on the region of Chamonix 
Mont-Blanc, France from 06 January 2009 to 24 March 2009 are available. A total of 8 SAR 
images are available every 11 days. Table [I] shows the main parameters of TerraSAR-X data. The 
area of interest covers the Argentiere glacier (Altitude: 2771m, 45.94628° N, 7.00456° E). The 
size of the domain is approximately 5km x 6km. Figure [4] shows the location and the image 
of Argentiere glacier captured on 06 January 2009. In order to obtain square pixels resolution, 
multi-look number of 5 for slant range and 3 for azimuth direction was applied. 
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Fig. 4. (Top) Location of the TerraSAR-X acquisition in the French Alps. (Bottom) The crop image on the Argentiere glacier 
area. The snowpack stratigraphic profiles calculated by Crocus snow model are given for 3 different altitudes on the Argentiere 
glacier: 2400m, 2700m and 3000m. The rectangles show the approximate positions of these altitudes on the TerraSAR-X images. 



For this study, meteorological forcing data provided by SAFRAN at 2400, 2700, and 3000 m 
altitude on horizontal terrain were used to drive the detailed snowpack model Crocus throughout 
the whole season 2008-2009 (starting on August 1st 2008). In order to carry out the comparison 
between the backscattering coefficients a S i m (obtained from executing the EBM using Crocus 
snowpack profile as input) and o~tsx (obtained from TerraSAR-X reflectivity), they need to be 
representative of the same area. Therefore we need to estimate the backscattering coefficients 
that well-represent the SAR reflectivities of the studied areas. The characteristic of a snowcover 
surface texture is spatially heterogeneous due to its strong variations of physical properties. A 
Gaussian distributed SAR texture hypothesis is therefore invalid. In recent studies, it has been 
proven that the texture of a SAR image of a heterogeneous medium can be modeled using the 



Fisher probability distribution [21 1, [22]. From the parameters of Fisher probability distribution, 
we can calculate the theoretical mean value which represents the backscattering coefficient of 
an area. In this study, the representative values of the backscattering coefficients of SAR image 
data on each altitude are obtained by calculating the mean value of Fisher-distributed texture of 
three regions (Fig. [4]) as in p2] [. 
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Frequency in Ghz Frequency in Ghz 



Fig. 5. Comparison of the contribution from surface, ground and volume backscattering mechanisms plotted in function of 
frequency /. The roughness parameters of surface and ground are: Left: gh = 0.4cm and I — 8.4cm corresponding to slightly 
rough and Right: — 1.12cm and I — 8.8cm corresponding to very rough (data taken from |23|). The volume backscattering 
coefficient is simulated using the snow profile in Table [TT] 

The roughness parameters of surface air-snow and ground are not available in the guess data 
calculated by Crocus, therefore the empirical values of the correlation length and the rms height 
have been taken from the measurements of Oh et. al. [23]. As we can see in figure |5} in the 
high frequency range (more than 10GHz), the contribution of surface and ground contribution 
are considerably low compared to the volume backscattering, regardless of slightly rough or 
very rough surface. Therefore we only concentrate on the volume contribution on the tests of 
data analysis method. According to the nature of a dry snow surface, the values of cr^ = 0.4cm 
and / = 8.4cm with a Gaussian type of surface spectrum, which correspond to a slightly rough 
surface, are used for the modeling in this study. With these surface and ground parameters set 
to constants, the original input vector x = [xcy OCMS x s x g ] in our case contains only the physical 
parameters of each layer of snowpack, which has the following form: 

X = [V-Crocus] = [Xl, X 2 , • ■ ■ , X^nf = [d\, d 2 , ■ ■ ■ , d n , pi, pi,..., Pnf (19) 

where di and pi are respectively the optical grain size and the density of i th layer of the snowpack. 
On the first iteration of the algorithm, x = x g and x 9 is given by the Crocus snow profile. 
The covariance matrix B, which represents the error of the input profile, i.e. of the Crocus 
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calculation, is a square (2n x 2n) definite positive matrix. Each element of B is computed as: 

Bjj = (7i.aj.jij (20) 



where a, = J E{(ei — Si) 2 } represents the standard deviation of error while calculating In our 
case, all element of x is estimated using the snow metamorphism model Crocus, therefore the 
variances of error are the same, which are experimentally estimated to 0.3 mm and 65 kg/m 3 
for the optical grain size calculation error and density calculation error respectively. 

The coefficient jij represents the correlation between errors of xi and xj and are modelled 
as: 

Jij = (3e- aAh *i (21) 

where Ah^ is the distance in cm between layer i and layer j. The values of a and (3 depend 
on different types of correlations and can be splitted into 3 cases: 

• Correlation d - d: a — 0.11 and (3 = 1 

• Correlation p - p: a = 0.13 and (3 = 1 

• Correlation d - p: a = 0.15 and (3 = 0.66 

These values are issued from an ensemble of slightly perturbated Crocus runs, i.e. obtained by 
differences in their meteorological inputs, over one winter season. The deviation between these 
runs, considered as elementary perturbations, have been then statistically studied and fitted with 



the eq. 21 model for the two considered variables and their crossed value. 

In this case study, the SAR data is available for HH channel, therefore the error covariance 
matrix R is a scalar which is equal to the variance of SAR image intensity on the studied area. 
The calculations of the variance on the three altitudes of Argentiere glacier gives the average 
value of R = 0.03. Nevertheless, after testing with different values, it has been observed that the 
output of the analysis algorithm is not very sensitive to this error factor. The scalar multiplication 
of 10 to 20 times the values of R doesn't show noticeable effect on the result. 



B. Results and Discussion 

Crocus snow stratigraphic profiles have been computed for 3 different altitudes over the 
Argentiere glacier, on the dates of the TerraSAR-X acquisitions. Two in-situ snowpack profiles 
measurements are also available at the altitude of 2700m on 30 January 2009 and 17 March 



November 15, 2012 



DRAFT 



16 



TABLE II 

Snow stratigraphic profile obtained from an in-situ measurement on Argentiere glacier on 30 January 2009 at the altitude of 

2700m. 



Snow depth 


Thickness 


Grain size 


Density 


(cm) 


(cm) 


(mm/10) 


(Kg/m 3 ) 











13 


13 


5 


210 


25 


12 


5 


290 


31 


6 


5 


310 


56 


25 


7.5 


220 


85 


29 


7.5 


300 


92 


7 


15 


340 


125 


33 


15 


430 


135 


10 


15 


370 


190 


55 


15 


430 



2009, and an example is shown in table |nj The level of liquid water content per volume at 
the time and location of measurement is below 2 percent, which means the snowpack can be 
considered as dry snow. Fig. [4] shows the approximate locations of each study area on the glacier. 

Fig. [6] shows the backscattering coefficients obtained over a period of time from three different 
methods: TerraSAR-X reflectivities, Crocus simulated profiles and simulation of analyzed Crocus 
profiles. Overall, the differences between the SAR reflectivities and the output of EBM simulation 
using Crocus initial guess profiles are approximately 2 to 6 dB. The EBM can overestimate the 
loss of EMW intensity while propagating through the snowpack medium due the assumption that 
snow particles are of spherical shape. The definition of snow optical grain size, the calculation 
of the effective permittivity and the phase matrix are also based on the same hypothesis. This 
assumption does not always hold in the natural environment where snow particles can have 
various shape and size. It is necessary to develop a more sophisticated method of modeling 
interaction between EMW and snow particles of different geometry properties. 
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Fig. 6. Results of simulation and analysis using a time series of TerraSAR-X acquisitions in 2009 and the corresponding 
Crocus output . The otsx (red) are the mean values obtained from Fisher probability distribution on a region of SAR image. 
The (Jsim (blue) represents the output of simulations using Crocus snowpack variables as input. The simulations again with the 
parameters after the data analysis process are shown in green. 



It can be observed that the gap between the TerraSAR-X backscattering coefficient and the 
simulation result of assimilated snow parameters is reduced to less than 1 dB. This shows that 
after the modification made by 3D-VAR, the analysed snowpack stratigraphic profiles give results 
of simulation closer to the backscattering value observed from radar. Some of the analysed o"o are 
less closed to the observation than the others (19 February 2009 of 2700m altitude or 2 March 
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Date: 13 Mar 2009. Altitude: 2700m. 



Fig. 7. The results of 3D-VAR data analysis method on some Crocus profiles. Each row contains 3 graphs that show the 
changes made by the analysis algorithm to each layer. Respectively from left to right: the snow optical grain size in mm, density 
in g/cm 3 and the backscattering coefficient of each layer in dB (output of the forward model) using initial Crocus profile and 
analyzed profile. From up to down are the graphs for the profiles on the altitude of 2700m, of which two in-situ measurements 
are also available: on 30 January 2009 (plotted on the same graph with Crocus profile on 28 January 2009) and 17 March 2009 
(same graph with Crocus profile on 13 March 2009). 



2009 of 3000m altitude). This may due to two reasons. First, the gradient descent using Newton 
method can converge to local minimum instead of global minimum. With a Crocus profile of 
50 snow layers, the algorithm involves of balancing 100 snow parameters in order to find the 
minimum of the cost function. The probability of having many local minima is significant high. 
Second, the modeling of the covariance matrix B needed to be further developed to have a more 
accuracy estimation of guess parameters' errors. Future works need to address these problematic 
using different kind of optimization method that can reduce the effect of local minimal and 
developing a better model for the error covariance matrix B. 
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TABLE III 

Comparison of bias and root-mean-square deviation (RMSD) between initial Crocus profiles and analysed profiles, with 

respect to the in-situ measurements 



Date 


Parameter 


Profile 


Bias 


RMSD 


28 Jan 2009 


Grain size (mm) 


Crocus 


0.43 


0.59 


Analysed 


0.45 


0.57 


Density (kg/m 3 ) 


Crocus 


110 


120 


Analysed 


40 


50 


13 Mar 2009 


Grain size (mm) 


Crocus 


0.49 


0.64 


Analysed 


0.46 


0.61 


Density (kg/m 3 ) 


Crocus 


120 


130 


Analysed 


50 


70 



Fig. [7] shows the detailed of the modifications of snow stratigraphic profiles done by data 
analysis process. The input parameters contain the snow optical grain size for Crocus, visually 
estimated grain size for the in-situ measurements, and the density of each snow layer. It can be 
observed that the modifications occur mostly on the near-surface layers. This can be due to two 
reasons: 

• The EMW at higher frequency has lower penetration rate. Depends on the compactness of 
the snowpack environment, X-band EMW can penetrate from 80 to 120 cm. This means 
the radar has little sensitivity to the characteristic of snowpack in the deeper layers. The 
EBM has taken into account this penetration rate through the calculation of attenuation. 
The deeper EMW penetrate, the higher value of attenuation is accumulated, and therefore 
the backscattering coefficient of the snow layers decreases exponentially from the surface 
layer to the ground layer. This can be observed from the graphs on third column of Fig. [7] 

• The error covariance matrix of measurements B is calculated based on the error correlation 



among layers. This correlation is based strongly on the distance between the layers (21). 
Large distance between two layers results in low value of correlation. Therefore the modifi- 
cations of the snow parameters of the near surface layers are considered independent from 
the deeper layers. 

Two in-situ measurements were carried out on 30 January 2009 (Tab. [TTJ) and 17 March 2009. 
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The stratigraphic profiles are plotted on the same graphs as the Crocus profiles of 28 January 2009 
and 13 March 2009 (Fig. [7]). The total simulated backscattering coefficients of these profiles are 
also displayed on Fig. |6j We calculate the bias and root-mean-square deviation (RMSD) between 
the initial (open-loop) Crocus profiles and the in-situ measurements, and compare to the bias and 



RMSD between the analysed profiles and in-situ measurements. Table III shows the comparison 
of these quantities. The results show the bias and RMSD between the analysed snow density 
and the measurements are much smaller than the initial guess (Crocus) snow density, hence the 
modifications made by the data analysis tend to approach the in-situ measurement. The analysis 
however shows little improvement with the modifications on the snow optical grain size, due to 
two reasons. First, the weight of the snow optical grain size in the covariance error matrix B, 
as well as in the adjoint model, is bigger than the weight of snow density. Hence it can also 
be noted that the values of optical grain size are not modified as much as the density (Fig. [7]). 
Second, the grain size used in Crocus is the snow optical radius, where the in-situ measurement 
uses the visually estimated grain size. Therefore the two quantities are not directly comparable 
to each other. 



V. Conclusion 

The results of this study show the potential of using data analysis method and the multilayer 
snowpack backscattering model based on the radiative transfer theory in order to improve the 
snowpack detailed simulation. The new backscattering model adapted to X-band and higher 
frequencies enables the calculation of EMW losses in each layer of the snowpack more accurately. 
Through the use of 3D-VAR data analysis based on the linear tangent and adjoint operator of the 
forward model, we have the possibility to modify and improve the snowpack profiles calculated 
by the detailed snowpack model Crocus. The output of this process shows that the discrepancies 
between the simulated profile and the in situ measurements are smaller after assimilation, and 
therefore could be further developed and used in real application such as snow cover area 
monitoring on massif scale or snowpack evolution through a period of time using series of 
spaceborne SAR image data. 

Future studies will be concentrated on developing the assimilation process. The 3D-VAR 
algorithm needs to be intergrated into Crocus, which means the analysed parameters of each 
step will be used as the input for the next step of initialization of Crocus. The result will be 
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an intermittent assimilation process where the snow stratigraphic profile generated by Crocus is 
continuously analysed and adjusted using TerraSAR-X data. 
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